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Abstract 



Tj- . We present a new numerical method for calculating the transfer of ionizing radiation, called 

C 2 -Ray (Conservative, Causal Ray-tracing method). The transfer of ionizing radiation in 
' diffuse gas presents a special challenge to most numerical methods which involve time- and 

spatial-differencing. Standard approaches to radiative transport require that grid cells must 
be small enough to be optically-thin while time steps are small enough that ionization fronts 
do not cross a cell in a single time step. This quickly becomes prohibitively expensive. We 
have developed an algorithm which overcomes these limitations and is, therefore, orders of 
magnitude more efficient. The method is explicitly photon-conserving, so the depletion of 
ionizing photons by bound-free opacity is guaranteed to equal the photoionizations these 
photons caused. As a result, grid cells can be large and very optically-thick without loss 
of accuracy. The method also uses an analytical relaxation solution for the ionization rate 
equations for each time step which can accommodate time steps which greatly exceed the 
characteristic ionization and ionization front crossing times. Together, these features make 
it possible to integrate the equation of transfer along a ray with many fewer cells and time 
steps than previous methods. For multi-dimensional calculations, the code utilizes short- 
characteristics ray tracing. The method scales as the product of the number of grid cells 
and the number of sources. C 2 -Ray is well-suited for coupling radiative transfer to gas and 
N-body dynamics methods, on both fixed and adaptive grids, without imposing additional 
limitations on the time step and grid spacing. We present several tests of the code involving 
propagation of ionization fronts in one and three dimensions, in both homogeneous and 
inhomogeneous density fields. We compare to analytical solutions for the ionization front 
position and velocity, some of which we derive here for the first time. As an illustration, 
we apply C 2 -Ray to simulate cosmic reionization in three dimensional inhomogeneous 
cosmological density field. We also apply it to the problem of I-front trapping in a dense 
clump, using both a fixed and an adaptive grid. 



Key words: Cosmology: theory; Galaxies: formation; galaxies: high-redshift; Intergalactic 



Preprint submitted to Elsevier Science 



5 February 2008 



medium; radiative transfer; methods: numerical 



1 Introduction 



The interplay between matter and radiation plays a crucial role in many astrophys- 
ical processes. The radiative feedback effects on hydrodynamic flows are partic- 
ularly strong for the case of hydrogen- and helium-ionizing radiation. The most 
dramatic manifestation of such radiative feedback, with far-reaching consequences 
for the present-day universe, was the reionization of hydrogen between redshifts 
z ~ 30 and z ~ 6 and helium by z ~ 3 — 4. The formation of the first ion- 
izing sources in the universe created expanding intergalactic H II regions. These 
eventually overlapped, leaving the universe largely ionized, as demonstrated by 
the absence of a Gunn-Peterson trough in the spectra of high-redshift QSO's and 
galaxies. On smaller scales, the change in pressure due to the heating associated 
with photoionization drives powerful flo w phenomena, such as photoevaporatiq n 



of minihalos during cosmic reionization (Shapi ro et al.l 12004: Il iev et al.ll20 05d). 
trigger ing and regulating star formation in Damped Lyman-a systems dlliev et al 
l2005ah . the formatio n of pillars in H II reg ions, evaporation of planetary disks 
in the Orion neb ula (Henney & Arthur 1998), the formation of planetary nebulae 



(Mellem alll994h. and other ne bulae, such as the famous ring around SN1987A 



(Chevalier & Dwarkadas 1995). 



Precise modeling of processes of radiative feedback is thus very important in or- 
der to understand all these phenomena. However, the development of effective and 
robust numerical implementations of radiative feedback, particularly in the case of 
optically-thick media, presents huge challenges. The full radiative transfer prob- 
lem adds to the usual three spatial and three velocity coordinates also angular and 
frequency dependencies, leading to a complicated, multi-dimensional problem. It 
also adds non-locality to the fluid equations, since distant sources can affect local 
dynamics, and introduces new requirements for the maximum size of the numerical 
time steps and cell sizes. Any naive, brute-force attempts to solve such problems are 
currently beyond the capabilities of even the fastest computers. The development 
of more efficient and robust radiative transfer algorithms is thus crucial. 

In the last few years there has been an intense development of numerical radiative 
transfer methods suitable for cosmology. However, compared to the much more 
sophisticated and mature methods developed for N-body and gas-dynamical sim- 
ulations, the modeling of radiative transfer effects is still in its infancy. Most of 
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the current state-of-the-art cosmological radiative transfer codes, with very few ex- 
ceptions, transport the radiation on pre-computed density fields, thus ignoring the 
dynamical back-reaction of the gas. While such approaches have their legitimate 
uses, the lack of gasdynamical feedback severely limits the questions that can be 
answered. Additional important limitations on all cosmological radiative transfer 
codes at present are imposed by their low resolution, which currently rarely ex- 
ceeds 128 3 computational cells in three dimensions. 



There are two basic classes of computatio nal radiative transfer methods currently i n 



use, moment methods ( Gnedin & Abel 



2001 



and ray-tracing metho ds (Mellem a et al 



1993 E 



1998 



Cenll2002l:lHaves & Norma5l2 003) 
Razoumov & Scottl 19991: 1 Abel et al 



Ciardi et al. 2000; Nakamoto et al. 2001; Sokasian et al. 2001; Razoumov et al. 



2002; Lim & Mellema 2003; Ma selli etalJ l2003: Sh apiro et alJ l2004: Bol ton et all 
2004: llliev et al.ll2005dl) . Each class of methods has its own advantages and disad- 
vantages. Generally, moment methods are fast and largely independent of the num- 
ber of ionizing sources, but are also fairly diffusive, which can l ead to incorrect 
results in some situations like e.g. prod ucing incorrect shadows (iGnedin & Abel 



2001 ; Cen 2002; Hayes & Norman 2003). On the other hand, ray-tracing approaches 



can be very accurate, but care should be taken to cover properly the space with 
rays, which often makes them computationally-expensive. This limits the number 
of sources that can be handled, and complicates the coupling of such methods to 
gas- and N-body dynamics. A wh ole new approach which may combine the best of 
both worlds, is being explored by Ri tzerveld et alJ (12003). 



The coupling of radiative transfer and hydrodynamics to model ionization fronts (I- 
fronts) is difficult because accurate finite-differencing in time and space generally 
requires time steps smaller than the cell-crossing time of the I-front (which can be 
highly supersonic). In addition, cell sizes must be small enough to be optically thin 
to ionizing radiation prior to the passage of the I-front. If cell sizes are limited to 
be optically thin, the Courant time drives the time step way down and the number 
of time steps way up. To make matters much worse, if the usual Courant condition, 
which limits the time step size to be less than the sound crossing time of a cell, has 
to be replaced by an "I-front Courant condition" involving the I-front speed rather 
than the speed of sound, this problem is greatly exacerbated. It is this combination 
of small cell size and small time step that makes the coupling of radiative transfer 
of ionizing radiation to gas and gravitational dynamics so computationally difficult 
by traditional finite-differencing approaches. In this paper we present a method that 
relaxes these requirements on both the spatial- and time-resolution. 



Previously we have developed an adaptive-grid (AMR), axisymmetric radiative- 
transfer and hydrodynamics code CORAL, which also contains non-eq uilibrium 
chem i stry and cooling due to H, He, and metal s (C, O, N, S and Ne) dRaga et al. 



1997: iMellema et al.lll998t IShapiro et alJl2004h . The latest versions of this code 
properly track fast, R-type I-fronts as well as slow, D-type I-fronts, but can require 
fairly large number of time steps in order to assure accuracy. We have applied this 
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code to a variety of astrophysical and cosmologic al problems, includin g photoe- 
vaporation of dense clumps in planetary nebulae (Mellemaetal. 1998) and pho- 
toevaporation of c osmological minihalos during reionization (Shapir o et al.ll2004l 



Iliev et al.N 2005d). The improved, three-dimensional (3D) successor of CORAL, 
different versions o f which use either AMR or uniform grid, is known as Yguazu 
(Rag a et al. lEoOOb). The uniform grid version of Yguazu was used e.g. to calculate 
pre-ionizatio n from a jet bow sh ock fRaga et al. 1999) and stellar jets moving into 
H II regions dRaga et al. 2000al) . iLim & Mellema (2003) used the to same basic 
methodology as in Yguazu, but in a different implementation, to study the interac- 
tion of two photoevaporating clouds in 3D. With adaptive mesh refinement these 
methods can reach effective resolution of 512 3 , or more. 



In order to improve upon the current implementation of radiative feedback in the 
CORAL and Yguazu codes, we present here a new approach to the photoionization 
calculations. Ensuring a high level of photon conservation helps relax the spatial 
resolution requirements of the code. All published ray-tracing methods have strong 
constraints on the time step in order to conserve photons. This makes combined hy- 
drodynamics and photoionization calculations expensive. Our method relaxes these 
constraints on the time step. The ultimate goal is to combine it with a hydrodynam- 
ics method, and hence speed and efficiency are essential. In the interests of length, 
in the current paper we will the describe our photoionization calculation and ray- 
tracing method, without discussing its coupling to hydrodynamics, which we will 
present in a follow-up paper. Our method is in fact also useful for 'stand-alone' or 
post-processing photoionization calculations, and that is how it is presented here. 



The structure of this paper is as follows. In § 2 we present our photon-conserving 
method for transferring the ionizing radiation and calculating the ionization rate. 
We also present a relaxation scheme to advance the non-equilibrium ionization rate 
equations across a finite time step, which is not limited by the ionization time. In 
order to use the method in a multidimensional setting, we need to cast rays from 
the sources. We describe our causal ray tracing scheme in § 3 and Appendix A. The 
treatment of multiple sources is discussed in § 4. In § 5 we present the tests we have 
performed to verify our method, in both ID and 3D. Finally, in § 6 we present the 
first illustrative applications of our method. 



2 Conserving photons 



Consider a continuum radiation field produced by an ionizing source with a spectral 
energy distribution of L u , traveling through a gas with a frequency-dependent op- 
tical depth t v . The flux of hydrogen-ionizing photons arriving at a distance r from 
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the source is given by 

1 7 L e- T "M 



where hu t h = 13.6 eV is the ionization threshold of hydrogen. The exact expression 
for the local ionization rate at a distance r from the ionizing source fo r hydrogen 
atoms with a cross section for ionizing photons a v is ( Osterbrocklll989l) 



1 [La e~ Tv ^ 

Wr) = ^y " \ v du. (2) 



"th 



The optical depth is defined, as usual, as 

t v = a v N m , (3) 

where iVni is the column density of neutral hydrogen. The expression in Eq. (2) is 
exact only at a given point in space and moment in time. However, in numerical 
simulations both space and time are necessarily discretized into finite-size cells 
and finite time steps. For finite cells the expression in Eq. (2) needs to be finite- 
differenced in a correct manner to ensure explicit photon conservation, which we 
discuss next. 



2. 1 Spatial discretization 



In a spatially discretized volume (a 'grid'), each spatial element does not have a 
single distance to the source, but spans a certain range Ar. Taking one ionization 
rate to be representative for this range is an approximation that is valid only if 
the grid cells are limited in size so that each cell is optically thin to the ionizing 
radiation. Since radiative transfer is computationally expensive, in general limiting 
the cell size in this way is prohibitive. As a result, the spatial discretization is often 
coarse, with very optically-thick cells. The effect of the approximation is that the 
number of photons absorbed by a 'grid cell' is no longer equal to the number of 
ionizations calculated for that cell. In other words, photons are not conserved, and 
ionizatio n fronts wi l l not t ravel at the correct speed. This problem was previously 
noted bv lAbel et all dl999h . who suggested that a better approach would be to force 



the ionization rate inside each halo cell to equal the absorption rate per cell used to 
attenuate the radiation in the transport algorithm. We shall adopt this approach and 
develop it further as follows. 

Consider a spherical shell of central radius r and width Ar, filled with neutral 
hydrogen of number density nni- Let N(r — Ar/2) be the number of ionizing 
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photons arriving at the shell per unit time, and N(r + Ar/2) the number of photons 
leaving the shell. The difference between these two numbers gives us the number 
of photons (per unit time) which were absorbed in the shell. These photons ionized 
a fraction of the nuiV she \\ hydrogen atoms in the shell, where V she \\ is the volume of 
the shell. The photoionization rate is then given by 

r _ ^~f)-^ + f) ; (4) 

W-HlKshell 



with 



47T 

Khell = ~7T 



ArV 
~2~) 



ArY 

TJ 



(5) 



Defining the optical depth from the source to r — Ar/2 as r„, and the optical depth 
between r — Ar/2 and r + Ar/2 (i.e. the optical depth of the cell) as Ar^, we can 
re-write Eq. (5) as 



oc 



—r 77 , (6) 

hv nniKheii 



Taking the limit of Ar„ <C 1 and Ar <C r (i.e. low optical depth per cell and also 
the distance from the source to the cell much larger than the size of the cell), one 
retrieves Eq. (2). Since transport is radial, this formula is also valid if we are only 
considering a sm all part of the shell. Equation (6) is equivalent to Eq. (15) from 



Abel etal.l (fl999) if we identify their V^n with V^heii- 



In Figure 1 we illustrate the difference between the local photoionization rate, 
riocai(^) from Eq. (2) evaluated at the entrance point of a cell and the photon- 
conserving photoionization rate, T, from Eq. (6) over the finite cell. We plot T 
for optical depths per cell Ar = 0.1, 1, 10, and 100, and ri oca i(r) (which is in- 
dependent of Ar) vs. the optical depth up to the cell, r. All rates are calculated 
in gray-opacity approximation (photoionization cross-section independent of fre- 
quency). We clearly see that, as expected, the local photoionization rate agrees well 
with the photon-conserving rates only for optically-thin cells (Ar < 0.1), in which 
case the two curves are barely distinguishable. However, the two rates start to dif- 
fer significantly even for moderately optically-thick cells (Ar ~ 1) and are com- 
pletely different for highly optically-thick cells (Ar > 10). In practice, the main 
differences when using photon-conserving finite-differenced rates would arise for 
moderately optically-thick cells for which the optical depth from the source r is 
low to moderate. For a cell which is either shielded or self-shielded (i.e. either r or 
Ar is large), the photoionization rate would be very low, and hence even incorrect 
photoionization rates would have no appreciable impact on the ionized fraction. 
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Fig. 1. Photoionization rate from Eq. (2) evaluated based on the optical depth, r, up to the 
cell (dotted line) compared to the photon-conserving photoionization rates from Eq. (6) for 
At = 0.1, 1, 10, and 100 (solid lines; top to bottom) vs. optical depth r from the source to 
the cell. We normalized the rates to 1 for r <C 1 and At <C 1, assuming r » Ar and gray 
opacities (i.e. a v = constant). 

Another way to think of the photon-conserving photoionization rate in Eq. (6) is to 
notice that for thin shell, Ar <C r, and gray opacity {p v =constant) we have 

r _ ri ocal (r) - r loca i(r + Ar) 
Ar 



This shows that the fractional error introduced by using Eq. (2) evaluated at the 
shell inner edge radius to give the ionization rate for all atoms in the shell, is given 
by [ri oca i(r) — ri oca i(r + Ar)]/r = Ar and hence grows as Ar grows. This again 
demonstrates that ri oca i evaluated at the inner cell boundary is a good approxima- 
tion to the finite-cell photon-conserving rate only for optically-thin cells. 



2.2 Temporal discretization 



When using the photoionization rate Y to integrate the ionization equations over a 
time step At, one normally assumes it to be constant during this time step. How- 
ever, since the ionizations and recombinations in the cell will change the density of 
neutral hydrogen, the optical depths r v and Atj, will also change. The usual solution 
is to take At which is small enough so the optical depths do not change appreciably 
over one step. For fast moving I-fronts this demands very short time steps. Various 
authors use different criteria in determinin g what the time st ep should be: a fraction 
of the photon trav el time over a cell A r / c dAbel et al.lll999h . several times the ion- 
ization time r _1 (Bo lton et al.l l2004). a small fraction of the timescale of change 
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of the H I fraction, nHi/(driHi/dt), (which sets the speed of the ionization front, 
Sha piro et al . 2004). Such small time steps make for expensive calculations if high 



accuracy is to be achieved. 

Since the above rate T is really a time-averaged rate (we only know what went in 
and what came out of the cell, not what happened in detail within the cell), a time- 
averaged value for Ar u can be expected to relax the constraint on the time step. 
We illustrate this with the following simplified model. Consider an infinitely-thin 
parcel of hydrogen gas of number density n illuminated by a time-varying photon 
flux F v (t) . Neglecting for the moment the effect of recombinations and collisional 
ionizations, the ionization rate equation can be written as 

d?/Hi r (9 . 



where yui is the neutral fraction of hydrogen. If we define a time-averaged pho- 
toionization rate, (T), as 

t+At 

(t> = Xt J nm', (9) 



then the solution to Eq. (8) is identical to that of 

—rr = - r h/Hi- (10) 
at 



This shows that, if only photoionizations contribute to the rate equation, it is correct 
to treat the ionization rate as a constant over the time step, no matter how large At 
is, as long as we take its value to be the time-average of T over At. If recombina- 
tions and collisional ionizations are taken into account, then this is not necessarily 
true since photoionizations that happen early in the time interval are more likely to 
be canceled by recombinations than those that happen later. However, as we shall 
show, using the time-averaged flux only introduces appreciable deviations from 
photon-conservation when the time step becomes comparable to the recombination 
time of the cell, and even in such cases the approximation holds fairly well. 

Solving for the ionization fractions is in itself not entirely trivial since we are deal- 
ing with stiff partial differential equations. Still considering only hydrogen, then 
the evolution of the ionized fraction x = 1 — ym can be written as 

dx 

— = (1 - x)(T + n e C H ) - xn e a H , (11) 



where n c is the electron density (itself dependent on the ionization fraction x) and 
Ch and an are, respectively, the (temperature-dependent) collisional ionization 
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and recombination coefficients for hydrogen. In this paper we assume the On- 
The-Spot (OTS) approximation (i.e. that the recombinations to the ground state 
are locally reabsorbed), in which case the recombination coefficient a# is equal to 
a B = 2.59 x l(T 13 (T/10 4 K)~ a7 cm 3 s~\ the Case B recombination coefficient 
for hydrogen at gas temperature T. 

If one takes T, n e , Ch and an to be constant, this equation has an analytical solution 

x(t) = x eq + Oo - x eq )e"' / ' i , (12) 



with 



£i = l/(r + n e C H + n c a H ) 



(13) 



and 



x 



cq 



r + n c C H 
r + u c Ch + n e an 



(14) 



Schmidt- Voigt & Koeppen (1987) suggested using this solution, and iterating for 
the electron density. The approach can be expanded to include other at omic species 
such a s helium and metals . This method was succ e ssfully used by e.g. Frank & Mellema 



1994jh lRaga etaD dl997h : lMellema et all dl998h : lshapiro et alJ d2004h : llliev et all 



(2005d), to name a few examples. In order not to have to repeatedly calculate the 



integrals over the spectrum, they are tabulated as function of the optical depth at 
u th , the ionization threshold of hydrogen, which, together with the spectral shape 
completely determines these integrals. We refer the reader to the above papers for 
further details on this method for solving the non-equilibrium chemistry equations. 

In the current context this approach for following the chemistry evolution equations 
is particularly appealing since it provides us with an analytical expression for the 
time-averaged ionization fraction 



(x) = X eq + (X - Xeq)(l - e A</<1 )^ 



(15) 



which can in turn be used to find the time-averaged optical depth of the cell (At), 
as follows: 

(Ar u ) = (1 - (x))n H a v Ar, (16) 



where (x) is given by Eq. (15). 
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Fig. 2. Flow chart of our method for a given computational cell. 

In order to remain self-consistent, the value of the (x) should then be used to modify 
T by changing both the optical depth and the neutral density: 



r L u e-^ l- e -(^> d ^ 

hv (riHl)^shell 



(17) 



"th 



which leads to an iterative process. Furthermore, since n e plays a role similar to 
T, we use its time-averaged value (n c ) when evaluating the ionization Eq. (11). 
The optical depth between the source and the cell edge r v is replaced by its time- 
averaged value, by adding, in causal order, all the time-averaged ( At„) of the cells 
lying between the source and the cell under consideration. This total optical depth 
(evaluated at the ionizing threshold of hydrogen, t = T„ =i , th ), is used to look up 
the values of the integral 

r>(r )= [—^-—dv, 08) 
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in pre-calculated tables. Then the photoionization rate in Eq. (17) is given by T = 
r(r ) - r'(r + Ar ), where Ar = Ar^. 



In practice, the iteration for finding the ionization state of a given cell proceeds as 
follows (see flowchart in Figure 2): 

(1) Initialize the mean ionization state to the initial values (given by the previous 
time step or the initial conditions). 

(2) Find the column density between the source and the cell 

(3) Iterate until convergence in the neutral fraction(s): 

- Calculate time-averaged column density of the cell [Eq. (16)]. 

- Calculate the photoionization rate T [Eq. (17)]. 

- Calculate the mean electron number density, (n e ), based on the current mean 
ionization state. 

- Calculate new and mean ionization states [Eqs. (1 1)-(15)]. 

- Check for convergence. 

To illustrate the ability of our method to obtain the correct result with a small num- 
ber of time steps, we consider the problem of formation and evolution of an H II 
region around a single source in an uniform, initially-neutral density distribution, 
which has an exact analytical solution (this is Test 1 in § 5 below, see there for 
details on the numerical parameters and the setup). We use a one-dimensional grid 
of 256 cells, gray opacities and uniform time steps. The optical depth per cell when 
it is fully-neutral is r ce n = 11.5. We solve this problem two ways, using the time- 
averaged photoionization rate in Eq. (17) with 100 time steps, and using the non- 
time-averaged rates in Eq. (6) evaluated at the beginning of each time step using 
10 2 , 10 3 , 10 4 and 10 5 uniform time steps. Results are shown in Figure 3. The time- 
averaged photoionization rates solution is essentially indistinguishable from the 
analytical result, with errors much smaller than 1%. To achieve a similar precision 
using the non-time-averaged photoionization rates we need up to ~ 10 5 time steps, 
or a factor of a 1000 more than in our method. The method with ~ 10 4 time steps 
and no time-averaging still gives an acceptable, but clearly inferior solution, which 
is off by a few per cent from the analytical result. Using any smaller number of time 
steps and no time-averaging leads to a completely incorrect evolution, whereby the 
I-front initially propagates much more slowly than it should have, and producing a 
much smaller H II region. It should be noted that ultimately even in those cases the 
correct Stromgren sphere is reached, but only at much later times. 

The reason why the method without time- averaging of the optical depth would 
need many more time steps is that it requires the time step, At, to be shorter (ide- 
ally much shorter) than the cell-crossing time, t cross = Ax ce \\/vi, (i.e. that the 
1-front does not cross more than one cell per time step). When the 1-front is very 
fast this condition demands very short time steps. In particular problems, where 
the I-fronts either propagate more slowly, or the fast-propagation p hase is short- 
lasting, utilizing such short time steps could be feasible (see e.g. IShapiro et al. 
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12 3 4 

t/t 

/ rec 

Fig. 3. Spherical I-front for point source in uniform gas: (a)(bottom panel) radius (in units 
of the Stromgren radius rg) vs. time (in units of recombination time t rec ). Time-averaged 
photon-conserving photoionization rates from Eq. (17) for 10 2 time steps (top line) and 
the non-time-averaged (but otherwise also photon-conserving) rates from Eq. (6) evalu- 
ated at the beginning of each time step, for 10 2 , 10 3 , 10 4 and 10 5 time steps (bottom to 
top lines). Note that the top curves (time-averaged and 10 5 step non-time-averaged) are 
almost indistinguishable, (top panel) Ratios of the numerical solutions over the analyti- 
cal one, r num /r ana i y t. Top lines, largely indistinguishable, are the time-averaged solution 
and non-time-averaged solution for 10 5 time steps, bottom line is the non-time-averaged 
solution for 10 4 time steps. 



2004), although still computationally-expensive. In more general situations, e.g. 
a cosmological density field, I-fronts often propagate very fast at least somewhere 
in the computational domain. In such problems using time-averaging of the optical 
depths is indispensable for this type of radiative transfer method. 



3 Ray tracing 



The method described in § 2 is one-dimensional, or in other words, along a ray. In 
order to use it in three spatial dimensions, we need to trace rays across the computa- 
tional domain. Different ray-tracing methods exist and could be combined with our 
method, as long as the rays are traversed in a causal order away from the source. 
Here we use a method called 'short characteristics', in which the rays for each point 
are constructed from previously calculated neighboring points, instead of casting 
inde pendent rays to each cell ('long characterist ics'). This is the approach used 
bv iRaga et al. dl999h and lLim & M ellemal d2003l) . with minor modifications. This 
method scales with the number of cells in the computational mesh. We describe our 
method in detail in Appendix A. 
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Fig. 4. Flow chart for the case of multiple sources. 
4 Treatment of Multiple Sources 



When multiple sources of ionizing radiation are present, we start by calculating 
the optical depths from each source outwards by casting rays to each cell and us- 
ing the method of short characteristics, as described in Appendix A. As we loop 
through the sources (in random order) we add the contribution of each source to the 
total photoionization rate, and we save the solution for the (time-averaged) neutral 
fractions. The latter implicitly communicates the presence of other sources to the 
current source. After having treated all sources we apply the total photoionization 
rate to the whole grid, this time without tracing any rays. This procedure is re- 
peated until the final neutral fraction satisfies our convergence criterion. Note that 
in the multiple source case we do not test for convergence when treating individual 
sources, but do so only at the end, when we apply the total photoionization rates 
from all sources. Figure 4 shows the above procedure in the form of a flow chart. 

Since ray-tracing is done for each individual source, the method roughly scales 
with the number of sources times the typical number of iterations needed for con- 
vergence. The latter number depends on the density field and the amount of overlap 
between the H II regions, so it is hard to give a general number. Experience shows 
it to be typically less than 10. 

Source numbers, positions, luminosities, starting times and lifetimes could be read 
in from a file or determined internally by the code. Currently all sources have the 
same spectrum, but within the framework of our method it is easy to implement 
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several types of sources (e.g. Pop. II and III starbursts, QSO's) by creating and 
using separate lookup tables for the photoionization rates for each source type. 



5 Tests 



5. 1 Expansion ofH II regions about a single source 



We first consider the simplest case of H II region expansion, namely around a sin- 
gle source of ionizing radiation in a fixed density field. Consider a source emitting 
iV 7 ionizing photons per unit time igniting in an initially-neutral medium consist- 
ing of pure hydrogen with number density n H . The source propagates an I- front 
into the neutral gas. The I-front propagation speed is determined by the balance at 
the 1-front of the incoming flux of ionizing photons and the flux of neutral atoms 
flowing into the 1-front and getting ionized. Assuming for simplicity that the 1-front 
is "sharp", i.e. the width of the neutral-ionized transition region is small compared 
to the 1-front radius, then this balance is expressed by the I-front jump condition 2 : 

n H ^r = F, (19) 
at 



where n H is the number density of hydrogen atoms and F is the flux of ionizing 
photons at the I-front. F is equal to the photon output of the source per unit time, 
iV 7 , minus the photons lost to recombinations in the already ionized volume, Vj = 

4vrrf/3: 



Anr 2 



iV 7 - 



\ 



n 2 H CaBdV 



Vi 



(20) 



/ 



where C = (n 2 ) / (n) 2 is the volume-averaged clumping factor of the gas (for a 
uniform gas C = 1), which increases the recombination rate of the ionized gas. 
Combining Eqs. (19) and (20), we obtain 



dVj_ 
d7 



iV 7 - 



n 



2 H Ca B dV 



(21) 



Vi 



2 It should be noted that strictly speaking the I-front jump condi tion in Eq. (19) shoul d 
be modified if the I-front is moving at relativistic speeds, see e.g. IShapiro et all yfX)5b). 
The numerical radiative methods should also take special care not to allow superluminal 
motions. However, in most astrophysical problems this does not become a significant issue 
and thus for simplicity in this paper we limit ourselves to considering non-relativistic I- 
fronts. 
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The width of an I-front is determined by the mean free path of the ionizing photons 
at the front. The assumption of "sharp" I-front is generally well-justified for soft 
ionizing spectra, but not necessarily for harder spectra, in which case the transition 
from neutral to ionized gas at the front could get fairly wide (see Sh apiro et al. 
2004, for further discussion and specific examples). 

We performed four tests of single-source 1-front propagation in several density 
fields relevant to astrophysics, which are listed in Table 1 and described in more 
detail below. In § 5.1.1 and § 5.1.2 we discuss the results in ID spherical symmetry, 
and 3D, respectively, with the temperature kept fixed at 10 4 K. For simplicity, we 
assume that the gas consists of pure hydrogen. While here we consider only the 
case of a pure hydrogen gas, our non-eq uilibrium chemistry solver can also easily 
accommodate multiple species (see e.g. IMellema et al.l ll998: Sh apiro et al.l l2004). 
We also assume a gray opacity. In this limit the photoionization rate from Eq. (17) 
becomes proportional to the number of ionizing photons produced by the ionizing 
source per unit time: 



In this simple setting there are exact analytical solutions for both the I-front position 
and velocity for all density fields we consider in our tests. To our knowledge, some 
of these solutions are derived here for the first time. 

In all tests the computational domain size was chosen so that it is only slightly 
larger than the final H II region size in order to keep the resolution roughly similar 
between the different tests. The source position for the ID, spherical symmetry tests 
is at the origin, r = 0, while for the 3D tests the source is in the center of the grid. 
We run each test four times in ID and four times in 3D, all possible combinations of 
two different spatial resolutions ("coarse" one with 16 cells, i.e. Ax cell = a; b0 x/16, 
and "fine" one with 128 cells in ID, as well as the corresponding 32 3 and 256 3 
cells in 3D), and two different choices of time steps, At coarse = t evo i/10 and 
Atfine = t e voi/100, where t evo i is the total time for which we follow the I-front. 
For our coarse-resolution runs we chose these very low spatial and time resolutions 
in order to demonstrate the conservative properties of our method. There are many 
problems where such low resolutions are expected, e.g. during reionization when 
we try to follow the initial evolution of multiple H II regions in a large simulation 
volume, so it is important to establish how reliable are the simulation results in such 
situations. Our "fine" resolution runs also employ relatively modest resolutions in 
both time and space, which enables us to easily run the ID and the 3D simulation 
results at equivalent resolutions in order to directly compare the results from the 
two. In all tests the parameters we picked are such that the individual cells start 
very optically-thick (see Table 1), so as to thoroughly test our code in such most 
demanding situations. 




(22) 
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Table 1 

Parameters of the test problems described in the text. The box sizes listed are for the 3D 
tests. For the corresponding ID tests we used a box size which was one-half of the val- 
ues listed. TQ yCe ii ym i n and TQ fie u^ max are the minimum and maximum optical depth per cell 
(when neutral) at the Lyman limit of hydrogen at the corresponding resolution (lr=low res- 
olution run, hr=high-resolution run). 

Test 12 3 4 

iV^s- 1 ) 10 54 10 51 10 51 10 54 

n m (cm- 3 ) 1.87 x 10~ 4 0.015 (^) 3.2 (^|^) 2 1.87 x 10~ 4 (1 + z)f 

C 5 1 1 5 

tcvoi(Myr) 500 15 1 500 

Xbox (cm) 5 x 10 24 1.4 x 10 22 6 x 10 21 7 x 10 25 (comov.) 



t rec (Myr) 130 8.2 (^) 0.04 130(1 + z)^ 

r s (kpc) 563 1.86 

t evoI /W 3.85 1.84 26 (9^) " 2 3.85(1 + z) 

T ,cell,minfr 184 91 34 50 

T0,cell,minMr 23 11 4 6 

T ,ceii, m axk 184 2200 3800 184 

To.ceH.mos.hr 23 2200 470 23 

We define the I-front position to be at the point where the gas is 50% neutral. Note 
that if the I-front width (which is typically 10-20 photon mean-free-paths, and is 
thus dependent on the type of ionizing spectrum and the spectrum hardening ahead 
of the I-front) is unresolved (the front is "sharp"), then the exact value of the neutral 
fraction used to define the I-front position is irrelevant. However, when the I-front 
structure is resolved (the front is "thick"), then its position (but not its velocity) 
depends on the neutral fraction at which we define it. Therefore, in such cases the 
I-front position obtained numerically could be offset from the position given by the 
analytical solutions, which all assume "sharp" I-fronts for simplicity. Inside a cell 
the I-front position is found by linear interpolation. The numerical value for the 
I-front velocity, v i, was obtained by simple finite-differencing of its position and 
the time, according to 

^i\ti,ave) "7 7 ' (23) 



-3 



3 
10 



where t i>ave = (ti + ij_i)/2. Whenever the time for the I-front to cross a cell is 
larger than the time step, the average speed of the front over a single time step is 
numerically poorly defined by Eq. (23). In that case, the speed was evaluated over 
a larger interval of time, comparable or greater than the cell-crossing time. 
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Test 1 involves the simplest situation, an I-front propagating in a uniform, con- 
stant gas density n H with a constant clumping factor C. There is a well-known 
analytical solution for the expansion of such an I-front, which is determined by 
two parameters, the Strom gren radius r s and the recombination time t rec (see e.g. 
iDopita & Sut herland l2003h . which are defined as 



rs 



3AL 



A^a B {T)Cn\ 



1/3 



(24) 



and 



t rec = [Ca B {T)n H ] 1 



(25) 



The analytical expressions for the 1-front position, 77 and velocity, v Is as a function 
of time are then 



r/ = r s [1 - exp(-t/t rcc 
r s exp(-t/t r 



,1/3 



Vi 



3t r 



[1 - exp(-t/t rec )] 



2/3 



(26) 
(27) 



i.e. the H II region reaches a finite radius, r s , and zero velocity at t — > 00 (in 
practice, after a few recombination times), at which point the recombinations in the 
ionized volume balance the new photons arriving from the source. In physical units 
we choose typical values for cosmological I-fronts propagating during reionization, 
with the gas density equal to the mean density of the universe at redshift z = 9 and 
a source with ionizing photon production rate N p h = 10 54 s _1 (Table 1). 

In Test 2 we study the propagation of an 1-front from a source in the center of a 
singular, decreasing density profile n H = n (r /r), where n = 0.015 cm -3 is the 
gas number density at the characteristic radius r = 5 kpc. This test is related to the 
problem of an I-front propagating outward fr om a source in the center of a galactic 
halo with a Navarro, Frenk & White profile dNavarro et al Jl997h (assuming the gas 
follows the dark matter density profile). For this density profile Eq. (21) reduces to 



dr rj 



(28) 



where we defined L = N y /(ATrn r ), K = n r CaB = r /t TeCj o, where t rcCi0 = 
(noCae) -1 is the recombination time at the characteristic density n . Eq. (28) has 
an analytical solution, which for initial condition r(0) = is given by 



ri(t) = rs \ 1 + Lambert W 



— cxp 



r t 



r s t 



- 1 



rcc,0 



(29) 
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where r$ = L/K is the Stromgren radius for this test and LambertW(x) is the 
solution of the algebraic equation y(x)e y ^ = x, which can be calculated e.g. using 
readily available public software. 

In Test 3 we follow the propagation of an I-front in a density profile n H = n (r /r) 2 , 
with a flat core of gas number density n and radius r . This density profile is 
steeper than the one we consider in Test 2, and t he H II region evolution is qual- 
itativel y different, as we show below (see also iFranco et al.l Il990l: IShapiro et al. 
2005b," for detailed discussion of I-front propagation in power-law density pro- 
files). This test, with the dimensional parameters we have chosen, n = 3.2 cm -3 
and r = 91.5 pc, (Table 1) resembles the problem of inside-out ionization of a 
dwarf galaxy formed at redshift z = 9 by an ionizing source at its center. In this 
case Eq. (21) reduces to 

£ = £ + (30) 
at rj 



where we defined 



AL 4 



L = v hm = - — - -n r Ca B , (31) 
47rn ro 3 



which physically is the terminal velocity v\ im of the 1-front for r — > oo, and K = 
n r 2 Ca B = rl/t mc ^ com . Assuming that iV 7 > Airr^nlCaB/S i.e. source is strong 
enough to ionize more than just the core, it is clear that vj > for all radii and 
thus the 1-front will never stop, eventually reaching the constant terminal velocity 
Vn m . Equation (30) for arbitrary values of the parameters K and L has a complex 
analytical solution for the initial condition r(0) > r . However, one particularly 
simple solution is obtained when L = 0, i.e. iV 7 = 167rrQnQ (70^/3, and r(0) = r , 
in which case the radius of the I-front is given by 

r 7 = r (l + 2t/t reCiCorc ) 1 / 2 . (32) 



This is the case we will use in our Test 3. In calculating the column densities for 
this test we use the weightings in Eq. (A. 5) with r = e > 0, where e< 1. 

Finally, our Test 4 is the same as Test 1, but for a cosmological I-front propagat- 
ing in a uniform-density medium with a density equal to the time-evolving mean 
background density of the universe 3 . The ionizing source is switched on at redshift 
z = 9. The I-front evolution of cosmological 1-front in an IGM w ith mean volume- 
averaged clumping factor C has an exact analytical solution ( Sha piro & Girouxl 

3 The cosmological parameters we use for the cosmological tests in this paper are H = 
TOkms^Mpc -1 , baryon density n b = 0.043 and n = 1 - S1 A = 0.27. 
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Fig. 5. Ionization front evolution for Test 1 (H II region expansion in an uniform gas) in ID 
spherical symmetry, using coarse grid (left panels) and fine grid (right panels). Shown are 
the analytical solution (solid, black), and the numerical solutions for At coarse (dotted, red) 
and Atfi ne (dashed, blue) for the I-front position (middle panels), ratio of the numerical to 
the exact analytical position (top panels) and the I-front velocity (bottom panels). 



19871) . given by 

y (x) = Xe x/x [xEi(2, X/x) - Ei{2, A)] . (33) 



Here y = [ri tC (t)/rs,i] 3 , ?7 jC = rj(l + z) is the comoving radius of the I-front at 
time if — ti + t, r s ,i = [3iV 7 / (Cas^c) 2 )] 1 ' 3 is the initial Stromgren radius, n H , c 
is the mean comoving density of hydrogen (defined at epoch of source turn-on, i.e. 
the scale factor is normalized a, = 1), 

U 

X = = tiCa B nH,c, (34) 

free,! 



is the ratio of recombination time at present to the age of the universe when the 
source turned on, and Ei(2, x) = ^-dt is the Exponential integral of second 
order. Note that the exponent in the analytical solution in Eq. (33), e xti ^, is gen- 
erally very large, which could easily lead to numerical overflow problems. There 
is a similar exponent, but with negative sign, in the Exponential integrals and it 
is therefore better to numerically evaluate them together, as e xti ^Ei(2, rather 
than separately. 



5.1.1 ID Ionization Fronts 

Figure 5 shows our results for Test 1 (constant density). The numerical results 
match the exact analytical solution quite well, even at very low spatial and tem- 
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Fig. 6. I-front evolution for Test 2 (H II region expansion in r 1 density profile) in ID 
spherical symmetry. Same notation as in Fig. 5. 



poral resolutions. The I-front position is never off by more than a fraction of a cell 
size at the corresponding spatial resolution. The only exception is the low temporal, 
but high spatial resolution case where the H II region size is initially underestimated 
by a few per cent (or ~ 1.25 cell-sizes). However, even in this case the agreement 
quickly improves and the final Stromgren radius is excellently reproduced. The nu- 
merical I-front velocity is also in excellent agreement with the analytical solution, 
regardless of the time- or spatial resolution, although the coarse spatial resolution 
runs tend to overestimate the I-front speed slightly when the I-front slows down 
below vj ps 100 km s" 1 . This could be expected, since the I-front cell-crossing 
time for this velocity is t cross ps 500 Myr, i.e. the 1-front has essentially stalled and 
remains inside a single cell, making the velocity estimate unreliable. This is con- 
firmed by our higher spacial resolution runs shown in Figure 5 (right) in which case 
the cell-crossing time for the same I-front velocity is 8 times shorter and hence the 
velocity is correctly reproduced down to few tens of km s _1 . 



In summary, for this test, the good temporal resolution is more important in achiev- 
ing good agreement with the analytical solution than is high spatial resolution, al- 
though the results from our numerical method are quite acceptable in all cases. 

The reason for this behavior is that for time steps of order of the recombination 
time, our method overestimates the number of recombinations. In this test prob- 
lem this happens when, within one time step, the H II region almost reaches its 
Stromgren sphere and the I-front essentially stalls. In this case, some cells near 
the front remain neutral for a large part of the time step. Using the time-averaged 
electron density in this case overestimates the number of recombinations. A better 
result is obtained if the time step, At is significantly shorter than the recombination 
time, t rec , while here At coarse = 50 Myr = 0.4t rec . 



20 



t [Myr] 



1 



1 

t [Myr] 



Fig. 7. I-front evolution for Test 3 (H II region expansion in r 2 density profile) in ID 
spherical symmetry. Same notation as in Fig. 5. 
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Fig. 8. I-front evolution for Test 4 (cosmological H II region expansion) in ID spherical 
symmetry. Same notation as in Fig. 5. 

Our results for Test 2 (1/r density) are shown in Figure 6. In this case the density 
profile is singular, and hence the evolution is dominated by recombinations from 
the start, rather than just at late times as was in Test 1 . The H II region reaches the 
correct Stromgren radius, with errors ~ 1% in all cases regardless of the resolution, 
but once again the early evolution of the H II region is much better followed (gen- 
erally within 1-2% or better) in the runs with higher temporal resolution. However, 
even for very poor time resolution the 1-front position is correct within ~ 10%, and 
to better than 4% after the first few time steps. The reason for this small initial dis- 
crepancy is that the recombination time in the inner, high-density cells is very short, 
shorter than (or at most of order) At coarse = 1.5 Myr, and hence the numerical I- 
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Fig. 9. Photon conservation for Tests 1-4 in ID. Cases plotted are: 16 cells, 10 time steps 
(red, dotted), 16 cells, 100 time steps (blue, short-dashed), 128 cells, 10 time steps (green, 
long-dashed), 128 cells, 100 time steps (black, solid). 

front initially propagates slightly slower than it should. Nonetheless, thereafter the 
I-front velocity is correctly reproduced in all cases up to the Stromgren sphere. 

Test 3 (1/r 2 density) proved most difficult for our numerical method to follow 
at low temporal resolutions (Figure 7). The H II region is smaller than the exact, 
analytical one by 5-10% percent (and up to 20% at the very beginning). Within 
one such relatively long time step (At coarse = 0.1 Myr) the I-front crosses the 
density profile core and advances significantly down the steep density gradient. The 
recombination time in the core is t rec , C ore = 0.04 Myr and At coarse = 2.5t reCiCore . 
Our time-averaging procedure slightly overestimates the mean optical depth over 
the time step in such conditions. Interestingly, the 1-front exact analytical velocity is 
nevertheless perfectly reproduced. Thus for coarse time steps the I-front propagates 
at the correct speed, but with slightly offset position. In the high temporal resolution 
cases the analytical I-front radius is matched to better than a few percent after the 
first few steps, regardless of the spatial resolution employed. 

Test 4 (expanding universe) resulted in excellent agreement between our numerical 
solution and the exact one, for all spatial and temporal resolutions (Figure 8). This 
was in fact expected, since in this case the recombination time is always longer than 
even our coarse time step, i.e. t rec /At coarse = 2.6/(1 + z)l , especially at later times 
as the background expands with the Hubble flow. We note that, as a consequence of 
this continuous decrease of the gas density, the Stromgren sphere is never reached 
and the 1-front velocity remains very fast at all times, always above 10 3 kms _1 , 
i.e. the 1-front is of fast, R-ty pe. This is a generic prope rty of global cosmological 



I-fronts, as was first shown by Shapiro & Giroux (1987) 



Finally, in Figure 9 we show the photon conservation numbers, which we define as 
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Fig. 10. I-front evolution for Test 1 (H II region expansion in an uniform gas) in 3D, using 
32 3 cells (left panels) and 256 3 cells (right panels), otherwise same notation as in Fig. 5. 
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Fig. 11. I-front evolution for Test 2 (H II region expansion in r density profile) in 3D. 
Same notation as in Fig. 10. 



the ratio of the total of all new ionizations and recombinations divided by the num- 
ber of photons provided by the ionizing source per time step (right), or integrated 
over time (left) for Tests 1-4 in ID spherical symmetry. We see that, even at the 
coarsest space- and time-discretizations, photons are generally conserved to within 
a few per cent. When either the spatial- or the time-resolution is not the coarse one 
the photon conservation of our method is generally better than 1%. 
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Fig. 12. I-front evolution for Test 3 (H II region expansion in r 
Same notation as in Fig. 10. 
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Fig. 13. Ionization front evolution for Test 4 (cosmological H II region expansion) in 3D. 
Same notation as in Fig. 10. 

5.1.2 3D Ionization Fronts 



For the 3D tests we use the same approach as in § 5.1.1, i.e. constant temperature 
at T = 10 4 K and gray opacities. The resolutions used are a coarse one at 32 3 
cells and a fine one at 256 3 cells. We placed the source in the center of the grid. 
Hence, the radial profiles are calculated at radial resolution of 128 and 16 cells, 
respectively, and are thus directly comparable to the ID results in § 5.1.1. Again, 
for each case we use two temporal resolutions, At coarse = i evol /10 and At fine = 
^evoi/100. Results are shown in Figs. 10-14. The radial cuts are along the positive 
x-axis. The results are largely identical to the corresponding ones in ID spherical 
symmetry. The numerical solution slightly overestimates the initial 1-front velocity 
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Fig. 14. Photon conservation for Tests 1-4 in 3D. Cases plotted are: 32 3 cells, 10 time steps 
(red, dotted), 32 3 cells, 100 time steps (blue, short-dashed), 256 3 cells, 10 time steps (green, 
long-dashed), 256 3 cells, 100 time steps (black, solid). 

for low temporal resolution cases in Tests 1 (uniform density) and 2 (1/r density). 
The final Stromgren spheres in the same Tests 1 and 2 have the correct sizes, to 
within few per cent, except for low spatial resolution runs of Test 2, where it is 
off by about 10%, or slightly more than 1 cell size. Once again the low temporal 
resolution runs in Test 3 (1/r 2 density) somewhat underestimate the H II region 
radius (although the match improves markedly at later times, to better than 5%). 
The I-front velocity in all cases is in good agreement with the exact result. Similarly 
to ID runs, Test 4 (uniform, expanding universe) exhibits the best agreement with 
the analytical solution, regardless of resolution. 

The level of photon conservation, either per time step or integrated over time in the 
3D tests is similar to the one we observed in the ID spherical runs, typically within 
few percent (Fig. 14). The conservation is worst for the low spatial resolution runs 
of Test 3 (1/r 2 density), largely due to a slight non- sphericity of the I- front (see 
Figure 15) caused by the steep gradient of the density profile. Nevertheless, even 
in these cases the photon conservation holds to within ~ 4%. The high spacial 
resolution runs conserve photons to better than a fraction of a percent in all cases. 



5.1.3 Testing Against Explicit Ionization Front Tracking for a Cosmological Den- 
sity Field 

As we have described, the method presented here calculates the evolution of the 
ionized fraction at each point in space by solving the ionization rate equations. As 
a way to check our implementation of this method for a single source surrounded by 
a static inhomogeneous density field, we have used a different, simpler approach 
which makes approximations similar to the ones used in obtaining the analytical 
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Fig. 15. I-front evolution. Shown is the 0.5 ionization level contour in time-sequence 
from t = 0.1t cvo i to t = t evo \, every 0.1t evo i for Tests 1-4, as labeled. All plots are for 
the case with coarsest-resolution in both space and time. On all occasions the I-front re- 
mains spherical at all times except for slight non-sphericity in Test 3 (1/r density). In all 
higher-resolution cases the I-front remains perfectly spherical at all times. 

solution in Eq. (27). The approximations are as follows. The I-front speed in any 
direction is given by the I-front jump condition in Eq. (19), using the flux on the 
ionized side and the density of neutral atoms on the neutral side of the front at that 
location. This flux is calculated according to Eq. (20) which assumes ionization 
equilibrium on the ionized side of the front. This allows us to attenuate the flux by 
integrating the recombination rate per unit area between the source and the position 
of the front along any direction, instead of calculating the optical depth. Along 
any given direction, the only thing that makes it necessary for us to solve these 
equations numerically, rather than analytically, is the fact that the gas density varies 
along the ray with distance from the source and so we must do a quadrature to 
get the integrated recombination rate along the ray to the position of the front. 
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The position of the front versus time along each ray (independent of every other 
ray) is determined by finite time-differencing the integration of the I-front velocity 
obtained from the jump condition as described above. Here we briefly describe our 
numerical implementation of I- front tracking. 

We use a set of isotropically distributed rays emanating from the source, and dis- 
cretize each of the rays radially into segments, so that segment i covers all radii 
r such that r\ < r < r i+ ±, and the center of the segment is given by rf +1 , 2 = 
(rf + rf +1 )/2. This is a version of "long-characteristics" approach to ray tracing, 
different from the short-characteristics employed by C 2 -Ray (see Appendix A). The 
density rii = n(r i+ i/ 2 ) at the center of the ray segment is determined by tri-linear 
interpolation from the eight nearest cell centers on the mesh. Along each ray, the 
I-front jump condition implies that 

dr = Njjryt) 
dt 47rr 2 n//(r)' 

where iV 7 (r) is the number of ionizing photons per unit time arriving at the I- front, 

r 

iV 7 (r) = iV 7i o - 4tt J a B n 2 H {r)r 2 dr, (36) 
o 

where iV 7i0 = iV 7 (r = 0). The discretization of Eq. (36) yields 

i-l 

Ni = iV 7)0 - Oi B n 2 H jAVi, (37) 

3 



where AVi = 4ir(rf +1 — rf)/3. Equation (35) is solved numerically to obtain the 
evolution of the 1-front along each ray. At any given time, the ionized fraction on 
the original grid is set to zero if the distance from the source is smaller than the 
radius of the I-front along that direction, and is set to one otherwise. 

We compare the results from the two methods as applied to the case of a single 
ionizing source in a cosmological density field (see § 6.2 for more details on the 
simulations). Since the explicit 1-front tracking assumes sharp I-fronts and fixed 
temperature, for closest direct comparison with the results of our method for the 
same test we assumed an ionizing source with a soft black-body spectrum with 
effective temperature T cff = 20, 000 K. The gas temperature is fixed to T = 10 4 K 
everywhere. In Fig. 16 we show the positions of the I-fronts obtained by the two 
methods at times t — 0.5 Myr and 1 Myr. The dotted blue lines show the I-front 
calculated with full radiative transfer method, while the solid black lines show the I- 
front position according to the explicit 1-front tracking. The thin red contours show 
the underlying density field. The H II regions assume a characteristic "butterfly" 
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Mpc VI >c 



Fig. 16. Full radiative transfer results from our method (thick dotted blue line) vs. explicit 
I-front tracking (thick solid black line) at times t = 0.5 Myr (left) and t = 1 Myr (right). 
Also shown are the density contours of the underlying density field (thin red lines). 

shape since the I-fronts, after they escape from the dense vicinity of the ionizing 
source, propagate faster in the voids and slower along the denser filaments. The 
results from the two methods show excellent agreement, thus verifying our ray- 
tracing procedure. 



6 Some simple illustrative applications 

6.1 Plane parallel ionization front on an adaptive mesh 

Our radiative transfer and ray-tracing method are directly applicable to any type 
of rectangular mesh, not just a uniform one, but also nested and adaptively-refined 
(AMR) meshes. We have currently implemented it for a plane-parallel I- front on an 
AMR mesh, in addition to the fixed, uniform mesh case. 

To illustrate the performance of our method on an adaptive mesh, we ran a simula- 
tion of a plane-parallel ionization front overrunning a dense spherical cloud. This 
problem was calculated both with a fixed mesh version, and with an adapt ive mesh. 
The A MR method used is the one implemented in the Yguazu code dRaga et al. 
2000b). It refines on a point by point basis, based on a given set of quantities, deter- 
mined by the user (these could include e.g. pressure, density, ionized fractions of 
various species, etc.) and employs refinement steps of a factor of two. For this par- 
ticular simulation we used a 3D box with sides of 8 x 10 18 cm, filled with uniform 
gas of 5 cm~ 3 and containing a ten times denser cloud of radius 10 18 cm positioned 
in the center of the computational volume. We used five refinement levels, and the 
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Time = 70 years Time = 270 years 




0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 

x (pc) x (pc) 

Fig. 17. Plane parallel ionization front crossing a high density region. Top row: calculated 
on a 5-level adaptive mesh; Bottom row: calculated on a fixed mesh. The snapshots show 
the density in the central xy plane (black=high density, white=low density), with the red 
contour showing the ionization front (x = 0.5). In the AMR simulation the dark points 
indicate the computational mesh points. 

maximum resolution (at finest refinement level) is 128 3 . We refined on the gradi- 
ents of the ionization fraction and density. For comparison we also ran the same 
problem on a fixed, uniform mesh of 128 3 cells. 

In Figure 17 we show the comparison between the AMR and the fixed mesh sim- 
ulation results. Our AMR implementation correctly identifies both the I-front and 
the clump and refines around them. The 1-front position and shape are identical 
at both times in the AMR and fixed-mesh results. This test demonstrates that our 
method interacts correctly with the adaptive mesh and properly tracks the I-front. 
For this particular problem the AMR simulation was completed approximately 10 
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Fig. 18. 1-front propagation around a single ionizing source in a cosmological density field 
at times t = 0.1 Myr (left) and t = 0.5 Myr (right). Image is a cut through the source 
y-plane and shows the color-coded gas density field (the darker regions are denser). The 
current position of the I-front is indicated by the white line. See Movies 1-3 for an animated 
view of the growth of these ionization fronts in three planes. 

times faster than the fixed mesh simulation, demonstrating the power of mesh adap- 
tivity to speed up codes without loss of accuracy. Further tests and full description 
of the AMR implementation of C 2 -Ray will be presented in a subsequent paper 
which will discuss the direct coupling of our radiative transfer method with hydro- 
dynamics (|Iliev et al 



6.2 Reionization of a cosmological density field 



As a further illustration and testing of our code, we have applied our method the 
simulation of reionization of a fixed cosmological density field at redshift z = 8.85. 
The gas distribution was obta ined from a PM+TVD N-body and gasdynamical sim- 
ulation (Shamro et al. 2005a) performed using the cosmological dynamics code de- 
scribed in R vu et alJ d 19931) . The simulation box size is 0.5 h~~ l Mpc, the resolution 
is 128 3 cells, 2 x 64 3 particles. 

We ran two simulations, the first one with a single ionizing source in the com- 
putational volume, the other with multiple sources. For the single- source run the 
ionizing source is positioned at the highest-density point of the grid (which for vi- 
sualization purposes is moved to the middle of the grid using the periodic boundary 
conditions), with luminosity L = 10 9 L Q and black-body spectrum of effective tem- 
perature T eff = 50, 000 K, as appropriate for a massive Pop II star, corresponding 
to an ionizing photon production rate of iV 7 = 6.822 x 10 52 s _1 . For the multiple- 
source run the ionizing sources are chosen to be the 6 most massive halos found in 
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t = 90000.0 years (z = 8.848) 



t = 200000. years (z = 8.847) 
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Fig. 19. Ionization of a cosmological gas density distribution by multiple sources at time 
t = 0.05 Myrs (top left) t = 0.07 Myrs (top right), t = 0.09 Myrs (bottom left), and 
t = 0.2 Myrs (bottom right). Image shows the color-coded gas density field, cut through 
the middle of the box (y-plane). The current position of the I-fronts is indicated by the 
white lines. Movies 4-6 show the growth of these ionization in three planes through the 
center of the volume. 

the box, with masses above 4 x 10 7 M<n, which mass a t that redshift corresponds to a 
halo virial temperature 10 4 K (Iliev & Shapir oll2001 ). below which temperature the 
gas cannot cool efficiently and form stars easily. The halos in the simulation box 
were found using a friends-of-friends halo finder with a linking length of 0.25. The 
ionizing photon production rate for each source is constant in time and is assigned 
assuming that each source lives t s = 3 Myr and has a constant mass-to-light ratio, 
emitting a total of f 1 = 250 ionizing photons per atom during its lifetime, as is 
appropriate for massive stars (Ilie v et al 1 l2005ch : 



(38) 
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Fig. 20. Ionized gas fraction for (a)(left) single-source and (b)(right) multiple ionizing 
sources on a cosmological density field: mass-weighted, x m , (solid) and volume-weighted, 
x v (dashed) (bottom panels) and their ratio (top panels). 

where M is the total halo mass and m p is the baryon mass. The source cell is 
assigned to the densest cell inside the 3 3 cells centered on the halo position. For 
simplicity all the sources are assumed to ignite at the same time. 

In the single-source case, the the H II region is initially confined to the dense region 
around the source. Once most of the gas there is ionized, the I-front quickly escapes 
into the lower-density voids, while the denser filaments temporarily trap it and be- 
come ionized only slowly, creating a characteristic "butterfly" shape (Fig. 18). Later 
on the H II region becomes somewhat more spherical and eventually most of the 
computational box becomes ionized at t ~ 1 Myr. 

When multiple ionizing sources are present in the computational volume (Fig. 19), 
the topology of the H II regions changes and becomes much more complex. Ini- 
tially, individual H II regions form around each source, but due to the significant 
source clustering at high redshift these small isolated H II regions quickly merge 
into a few larger ones, with some remaining neutral pockets of denser gas along 
the filaments inside them which temporarily self- shield and trap the propagating 
I-front. The shape of the H II regions is quite non-spherical at all times, dictated by 
the interplay between the source clustering and the underlying density distribution. 
The computational domain becomes completely ionized within ~ 1 Myr. 

This shows that reionization is not as simple as either "inside-out"- dense regions 
ionized first, low-density regions last - or "outside- in" - low-density reg i ons, o r 
voids, ionized first, dense regions last, claimed by Miralda-Escude et al. (2000). 



The dense regions in the immediate vicinity of the sources were ionized first, but 
then the I- fronts quickly escaped into the voids, and thereafter the overdense and 
under-dense regions were being ionized simultaneously. This point is further il- 
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lustrated in Figure 20, which shows the mass- and volume- weighted ionized frac- 
tions for both the single- and the multiple- source runs. Initially, in both cases the 
mass-weighted ionized fraction is significantly larger than the volume-weighted 
one. However, the ratio quickly, within 0.1-0.2 Myr, drops to ~ 1. The evolution 
is somewhat different thereafter. In the single- source case the ratio becomes some- 
what less than one, as the voids become ionized faster than the dense regions. In 
the multiple-source case, on the other hand, the ratio stays very close to 1, indi- 
cating that while the voids (which take most of the volume) are ionized relatively 
quickly, the dense knots and filaments (which contain most of the gas mass) are 
also being ionized at the same time, which keeps the ratio of the mass- to volume- 
ioni zed fraction similar. W ith this simple example we thus confirm the results of 
e.g. lRazoumov et all d2002h . 



7 Conclusions 



We presented a new method, called C 2 -Ray, for photoionization calculations and 
radiative transfer in optically-thick media. Our method is explicitly photon-cons- 
erving, which allows us to relax the often-required condition for the computational 
cells to be optically-thin, and thus be able to obtain accurate results even for very 
coarse spatial discretizations. Furthermore, we introduced time-averaged optical 
depths, which we obtain from a relaxation solution of the non-equilibrium chem- 
istry equations, which allow us to use similarly coarse time-discretizations without 
any significant loss of accuracy. These developments result in a very fast, efficient 
and accurate calculation of the evolution of the ionization state of a gas distribution 
regardless of the spatial- and time-resolution. 

We tested our method systematically, and in significant detail, in several astro- 
physically relevant situations in both ID (in order to test its basic functionality) 
and in full 3D, at various resolutions in time and space. We compared the results 
against the corresponding exact analytical solutions, some of which we derived 
here for the first time. The results of these tests show that our method follows the 
analytical solutions to within a few percent, and conserves photons with the same 
accuracy or better, even when using extremely coarse grids and very long time 
steps. We have also shown that our radiative transfer method can be combined with 
adaptive mesh refinement (AMR), which leads to substantial increase in the cal- 
culation speed and much lower memory requirements. These tests demonstrate the 
importance of extensive testing of RT codes in a wide range of density fields - both 
for code validation and for possible adjustments of free parameters for better han- 
dling of special cases, e.g. steep density gradients. A project to thoroughly compare 
different codes used for cosmological radiative transfer methods is underway (Iliev 
et al. 2005, in preparation). 

Our radiative transfer method imposes no limitations on the computational cell size 
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and optical depth, and imposes no hard limits on the time step sizes. The only re- 
quirements we found are related to the desired accuracy of the solution obtained, 
for high accuracy the time step should be significantly shorter than the local re- 
combination time. Longer time steps, comparable to the recombination time, lead 
to gradual loss of accuracy, but no incorrect or unstable solutions. 



As a first simple application, we studied the reionization of a cosmological den- 
sity field by both a si ngle- and multiple- sources confirmed the findings of e.g. 



Razoumov et al. (2002) that neither inside-out, nor outside-in models describe the 



progress of cosmological I-fronts correctly. We showed that, instead, the mass- and 
volume-weighted ionized gas fractions are generally very similar, indicating that 
the dense and underdense regions are ionized at the same time. 



The accurate I-front tracking over long time steps and coarse spatial resolutions 
makes our code ideal for direct dynamical coupling to multi-dimensional gasdy- 
namics and N-body codes, where the evolution time-scales are generally much 
longer than the ionization and 1-front crossing times. Allowing for larger, optically- 
thick cells further relaxes the time step limits imposed on the gasdynamic evolu- 
tion through the Courant condition. In a follow-up paper we will address how we 
combine our C 2 -Ray method with 3D gas- and N-body dynamics, including more 
details on the coupling to AMR and calculating the effects of photo-ionization heat- 
ing. 



The C 2 -Ray's high computational efficiency and correct 1-front evolution over long 
time steps make it also a very attractive tool for studying reionization using high- 
resolution precomputed density fields. As a ray-tracing method, its performance 
scales linearly with the number of ionizing sources and the number of computa- 
tional cells, but our tests show that on the available hardware we can perform such 
post-processing calculations on large computational meshes of 400 3 to 800 3 cells 
and for thousands of ionizing sources. In a future paper we will describe the results 
of these calculations and show how we use them to make detailed observational 
predictions e.g. for the large-scale 21 -cm signal from the Epoch of Reionization. 
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A Ray tracing 



Here we describe how the column densities for a given cell are constructed using the 
short-characteristics method. For this we consider a point d with position (x, y, z) 
at mesh position (i, j, k), see Fig. A.l. The source point s in physical coordinates is 
(x s , y s , z s ) or in mesh coordinates (i s , j s , k s ). The point where the ray enters the cell 
is called c. Assuming that the cells are cubical, the ray calculation is most easily 
done in mesh coordinates, which is what we will use. 

For our method we need two column densities for each cell, the column density to 
the cell, called N c and the column density over the cell, denoted AiV. N c corre- 
sponds to the optical depth t u , and AN corresponds to Ar u in Eq. (6). 

AiV is simply given by the n^ds, with ds the path length ds through the cell, 
which is twice the distance between c and d and can be shown to be 



ds 



(i ~h) 2 + (] ~ 7s) 2 



in units of cell size. The appropriate value for V^heii m Eq. (6) is then given by 
4-7rr 2 d ds, where r sd is the distance between the source and the point d. 

The column density to the cell N c is constructed from the column densities of the 
neighboring cells. In order to pick the correct neighbors, we need to establish where 
the ray enters the cell around d. This can be done using the differences Ai = i — i s , 
Aj = j — j s , Ak = k — k s which show whether the ray from the source enters the 
cell through the x, y, or z plane. E.g. if \Ak\ > \Aj\ and \Ak\ > \Ai\ the crossing 
is on the constant z -plane. Considering this case we calculate the point c where 
the ray enters the cell, with mesh coordinates (i c ,j c , k c ), where k c = k — |crfc and 

lAfel 
a k ~ -Ak- 

The neighbors closer to the source are the four cells with mesh positions el: (i, j, k— 
a k ), e2: - aj, k - a k ), e3: (i - a^j, k - a k ), and e4: (i - ai,j - aj, k - a k ), 
with a i: j = ^pi. From these we construct the column density at the crossing point 
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c: 

N c = Wl N el + w 2 N e2 + w 3 N e3 + w 4 N e4 , (A.2) 

with Wi a set of normalized weights. These weights can be constructed in different 
ways. The most straightforward one is taking the x and y distances from point c to 
the corner points: 

5 i = 2|» c -t + o- < /2| 
^ = 2|j c -j + a,/2| 

and using 

^ = (1-^(1-5;) 

V 3! (A.4) 

w 3 = (1 - 5i)5j 



W4 = SiS 



.1 



With this weighting N c = N el if the ray is parallel to the z-axis, and iV c = N eA if 
the ray travels diagonally over the mesh. 

However, we found that this simple geometric weighting in case of a ver y clumpy 



medium gives too much spreading of shadows. As discussed previously in Rag a et al 



(1999), we suppress this by adding an extra optical-depth-dependent factor in the 



weights wi — w 4 . 

W 2 



majX.(T0,Te2)^2 n W n c\ 

(l-5i)6i ' 

w s = — r — t6 

W4 = -. 

max(ro,T e 4) 2^ n w n 



The value of tq was determined empirically at 0.6. For this number we get the best 
photon-conservation for general clumpy media, when there are no strong density 
gradients present. In all tests and applications presented in this paper we use the 
weightings in Eq. (A. 5), unless otherwise noted. 

In the special case of density field with strong density gradients (e.g. n oc r~ 2 , 
as in one of the tests we discuss in § 5) combined with low spatial resolution, the 
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weightings in Eq. (A. 5) with r = 0.6 do not produce good photon conservation. 
This is easy to understand since steep density gradients combined with very coarse 
resolution lead to incorrect representation of the density field on the grid and of 
significant grid-induced anisotropics, which are only worsened by the weightings 
in Eq. (A.5). In such cases it is better to use the weightings in Eq. (A. 5) with a 
small value for r (r = e > 0), i.e. essentially the weightings in this case are simply 
inversely proportional to the corresponding optical depths (or, equivalently, column 
densities), since max(r ,r e j) = r ej j for i = 1,4. When the spatial resolution is 
sufficiently high (e.g. 128 3 or 256 3 mesh in the tests discussed in § 5) both types 
of weightings work well and produce almost identical results, with very high level 
of photon conservation, generally better than 1%, regardless of the time resolution 
employed. 

For |Aj| > |A/c| and |Aj| > |Ai| we are dealing with a y-plane crossing and for 
\Ai\ > \Ak\ and |Ai| > |Aj| with an a; -plane crossing. For these two cases the 
calculations are done in a similar way as above, but with the y and x coordinates 
taking the role of the z coordinate, respectively. 

Both the short characteristic method (which uses the column densities from the 
neighbors, as described above), and the time-averaged optical depth calculation, 
require the ray-tracing to be causa l. To achieve this t he mesh has to be traversed in 



a particular order, as described in Raga et al. ( 1999): we first trace away from the 



source along the x-axis (y = y s , z = z s ), next we do strips of constant y parallel 
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to the x-axis. After the source plane (z = z s ) has been traced, we move to z -planes 
away from the source, and trace each plane in the same way as the source plane. 
Since each cell has its own ray segment, the short characteristics method scales 
with the number of cells in the computational mesh. 

A.0. 1 Shadowing by a Dense Gas Cloud and I-front trapping 

Some of the main effects due to radiative transfer in realistic simulations are slow- 
ing and possible trapping of the I-front in denser regions and the corresponding 
casting of shadows behind these dense regions. In this section we are testing our 
numerical scheme in such a situation. We use the initial conditions from Test 1 in 
Table 1, except that the computational box is 1 Mpc in size. We position a dense 
(gas number density n = 10 3 cm -3 ), uniform, rectangular slab with size 2 x 20 23 
by 5 x 10 23 by 5 x 10 23 cm at a distance 0.083 Mpc from the source (sizes are 
arbitrary and have no particular physical significance). The temperature is fixed at 
t = 10 4 K everywhere. Results for a 128 3 mesh are shown in Figure A.2. As ex- 
pected, the I-front is trapped inside the dense region within one cell. The shadows 
cast are clean and sharp, with only modest spurious diffusion in the shadow behind 
the clump. The I-front outside the shadow is not affected by its existence and keeps 
its original, perfectly circular shape. 
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Fig. A.2. Shadows test. Same parameters as in Test 1 above, but introducing a dense 
(n = 10 3 cm -3 ) uniform slab (red) for resolution of 128 3 cells. Source is in the mid- 
dle of the computational box (blue circle). Contours are of the 50% ionized fraction in a 
time-sequence from t = 20 Myr to 100 Myr, with intervals of 20 Myrs. The thicker line 
is for t = 20 Myrs and its thickness indicates the width of the ionization front. The three 
panels for each case are cuts through the center along x-y, y-z and x-z plane, as labeled. 
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